Solving the Bethe-Salpeter equation for a pseudoscalar meson in 

Minkowski space 



V. Sauli 

CFTP and Departamento de Fisica, Instituto Superior Tecnico, 
Av. Rovisco Pais, 1049-001 Lisbon, Portugal, and 

Department of Theoretical Physics, Nuclear Physics Institute, 
Rez near Prague, CZ-25068, Czech Republic 

Abstract 

A new method of solution of the Bethe-Salpeter equation for a pseudoscalar quark-antiquark 
bound state is proposed. With the help of an integral representation, the results are directly 
obtained in Minkowski space. Dressing of Green's functions is naturally taken into account, thus 
providing the possible inclusion of a running coupling constant as well as quark propagators. First 
numerical results are presented for a simplified ladder approximation. 
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I. INTRODUCTION 



Among the various approaches used in meson physics, the formahsm of Bethe-Salpeter 
and Dyson-Schwinger equations (DSEs) plays a traditional and indispensable role. The 
Bethe-Salpeter equation (BSE) provides a field-theoretical starting point to describe hadrons 
as relativistic bound states of quarks and/or antiquarks. For instance, the DSE and BSE 
framework has been widely used in order to obtain nonperturbative information about the 
spectra and decays of the whole lightest pseudoscalar nonet, with an emphasis on the QCD 
pseudo-Goldstone boson — the pion Moreover, the formalism satisfactorily provides a 
window to the 'next-scale' meson sector, too, including vector, scalar 0] and excited mesons. 
Finally, electromagnetic form factors of mesons have been calculated with this approach for 
space-like momenta 0. 

When dealing with bound states composed of light quarks, then it is unavoidable to 
use the full covariant BSE framework. Nonperturbative knowledge of the Green's function, 
which makes part of the BSE kernel, is required. Very often, the problem is solved in 
Euclidean space, where it is more tractable, as there are no Green's function singularities 
there. The physical amplitudes can be then obtained by continuation to Minkowski space. 
Note that the extraction of mass spectra is already a complicated task [3], not to speak of 
an analytic continuation of Euclidean-space form factors. 

When dealing with heavy quarkonia or mixed heavy mesons like Be (found at Fermilab 
by the CDF Collaboration [5|), some simplifying approximations are possible. Different 
approaches have been developed to reduce the computational complexity of the full four- 
dimensional (4D) BSE. The so-called instantaneous |6| and quasi-potential approximations 
0] can reduce the 4D BSE to a 3D equation in a Lorentz- covariant manner. In practice, such 
3D equations are much more tractable, since their resolution is less involved, especially if one 
exploits the considerable freedom in performing the 3D reduction. Also note that, contrary 
to the BSE in the ladder approximation, these equations reduce to the Schrodinger equation 
of nonrelativistic Heavy- Meson Effective Theory and nonrelativistic QCD [8|. However, the 
interaction kernels of the reduced equations often correspond to input based on economical 
phenomenological models, and the connection to the underlying theory (QCD) is less clear 
(if not abandoned from the onset). 

In the present paper, we extend the method of solving the full 4D BSE, originally devel- 
oped for pure scalar theories 0, to theories with nontrivial spin degrees of freedom. 



Under a certain assumption on the functional form of Green's functions, we develop a method 
of solving the BSE directly in Minkowski space, in its original manifestly Lorentz-covariant 
4D form. In order to make our paper as self-contained as possible, we shall next supply 
some basic facts about the BSE approach to relativistic mesonic bound states. 

The crucial step to derive the homogeneous BSE for bound states is the assumption that 
the bound state reflects itself in a pole of the four-point Green's function for on-shell total 
momentum P, with = M|, viz. 

G(4)(p,p', P) = V -'^^(p^Pos)^Ap'^Pos) J 



where Ep^ = yP^ + Mj and Mj is the (positive) mass of the bound state characterized by 
the BS wave function ijjj carrying the set of quantum numbers j. Then the BSE can be 
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conventionally written in momentum space like 



S{\p+, P)^Pip, P)S^\p^, P) =-ij ^V{p, k, P)^{p, P) , (1.2) 

p+ = P + aP , 

p„ = p — (1 — a)P , 



or, equivalently, in terms of BS vertex function F as 

r(p, P) = -^J P)Si{K, P)T{p, P)S2{k., p) , (1.3) 

where we suppress all Dirac, flavor and Lorentz indices, and a G (0,1). The function V 
represents the two-body- irreducible interaction kernel, and Si {i = 1, 2) are the dressed 
propagators of the constituents. The free propagators read 



) = /+!"^ . . (1.4) 



Concerning solutions to the BSE (11. 3p for pseudoscalar mesons, they have the generic 
form 12 

r(g, p) = 75[rA + Tbq.p 4+'^c r + '^D 4 r + '^E F (i.s) 

where the Fj, with i = A, B,C, D, E, are scalar functions of their arguments P,q. If the 
bound state has a well-defined charge parity, say C = 1, then these functions are even in 
q.P, and furthermore = — F^;. 

As was already discussed in Ref. 13|, the dominant contribution to the BSE vertex 
function for pseudoscalar mesons comes from the first term in Eq. (11. 5p . This is already true, 
at a 15% accuracy level, for the light pseudoscalars vr, K, r], while in the case of ground-state 
heavy pseudoscalars, like the r]c and rji,, the contributions from the other tensor components 
in Eq. (II. 5p are even more negligible. Hence, at this stage of our Minkowski calculation, we 
also approximate our solution by taking F = 75F^. 

The interaction kernel is approximated by the dressed gluon propagator, with the inter- 
action gluon-quark-antiquark vertices taken in their bare forms. Thus, we may write 

V{p, g, P) = g\K)D^,{p - g, K)Y ® Y , (1-6) 

where the full gluon propagator is renormalized at a scale n. The effective running strong 
coupling as is then related to g through the equations 

9\k)D,^{1,k) =«.(U)^-e/(K)-^, (1.7) 



l-n(g2,K) 



From the class of ^-linear covariant gauges, the Landau gauge ^ = will be employed 
throughout the present paper. 
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In the next section, we shall derive the solution for the dressed-ladder approximation 
to the BSE, i.e., all propagators are considered dressed ones, and no crossed diagrams are 
taken into account. The BSE for quark-antiquark states has many times been treated in 
Euclidean space, even beyond the ladder approximation. Most notably, the importance of 
dressing the proper vertices in the light-quark sector was already stressed in Ref. [3], so 
our approximations are certainly expected to have a limited validity. Going beyond the 
rainbow (7^) approximation is straightforward but rather involved. (For comparison, see 
the Minkowski study of Schwinger-Dyson equations published in Refs. fisl . [l6|). the latter 
paper including the minimal-gauge covariant vertex instead of the bare one). In the present 
paper, we prefer to describe the computational method rather than carrying out a BSE 
study with the most sophisticated kernel known in the literature. 

The set-up of this paper is as follows. In Sec. 2 we describe the method of solving the BSE. 
As a demonstration, numerical results are presented in Sec. 3. Conclusions are drawn in 
Sec. 4. The detailed derivations of the integral equation, that we actually solved numerically, 
are presented in the Appendices. 



II. INTEGRAL REPRESENTATION AND SOLUTION OF THE BSE 

In this section we describe our method of solving the BSE in Minkowski space. It basically 
assumes that the various Green's functions appearing in the interaction kernel can be written 
as weighted integrals over the various spectral functions (i.e., the real distribution) p. 

More explicitly stated, the full quark and gluon propagators, the latter ones in the Landau 
gauge, are assumed to satisfy the standard Lehmann representation, which reads 

s{i) = r dj-^^^tii±j!M ^ (2.1) 

JO V — uj + ie 



GAl) = dco n . P^uiO , (2.2) 



P-to + ie 

where /? is a real istribution. Until now, with certain limitations, the integral representations 
(12. ip and (12.21) have been used for the nonperturbative evaluation of Green's functions in 



various models [17j. However, we should note here that the true analytic structure of QCD 



Green's functions is not reliably known (also see Refs. [18|, ll9|, l20| ) , which studies suggests the 
tructure given by (12.11) and (12.21) is not sufficient if not excluded. In this case, the lehmann 
representation or perhaps the ussage of real p in the integral representation (12. ip and 
(12. 2p can be regarded as an analyticized approximation of the true quark propagator. The 
complexification of rho within the complex integration path is one of the straightforward 
and questionable generalization [2l[. The general question of the existence of Lehamnn 
represintation in QCD is beyond the scope of presented paper and we do not discussed the 
problem furthermore. 

Furthermore, we generalize here the idea of the Perturbation Theory Integral Represen- 
tation (PTIR) [9|, specifically for our case. The PTIR represents a unique integral repre- 
sentation (IR) for an n-point Green's function defined by an n-leg Feynman integral. The 
generalized PTIR formula for the ra-point function in a theory involving fields with arbitrary 
spin is exactly the same as in the original scalar theory considered in Ref. j^, but the spec- 
tral function now acquires a nontrivial tensor structure. Let us denote such a generalized 
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weight function by p{a,Xi). Then, it can be clearly decomposed into the sum 



pi^Oi, Xj)scalar theory ^ ^ ^ Pji,^i •^i)'^^ 



(2.3) 



j 



where a,Xi represent the set of spectral variables, and j runs over all possible independent 
combinations of Lorentz tensors and Dirac matrices Pj. The function pj{a, xt) just represents 
the PTIR weight function of the j-th form factor (the scalar function by definition), since it 
can obviously be written as a suitable scalar Feynman integral. Leaving aside the question of 
(renormalization) scheme dependence, we refer the reader to the textbook by Nakanishi [9| 
for a detailed derivation of the PTIR. The simplest examples of such "generalized" integral 
representations corresponds with Lehmann representations for spin half (12. ip and spin one 
propagators (12.21) . 

Let us now apply our idea to the pseudoscalar bound-state vertex function keeping in 
mind that the singularity structure (given by the denominators) of the r.h.s. of the BSE 
is the same as in the scalar models studied in Refs. [lo|, the appropriate IR for the 
pseudoscalar bound- state vertex function r^(g, P) should read 



where we have introduced a useful abbreviation for the denominator of the IR (12. 4p . viz. 



with a free integer parameter. 

Substituting the IRs ([23]), (j^II]) into the r.h.s. of the BSE (O), one can analytically 

integrate over the loop momenta. Assuming the uniqueness theorem we should arrive 
at the same IR (12.40 . because of the r.h.s. of the BSE (11.30 . The derivation is given in 
Appendix A for the cases A^ = 1, 2. 

In other words, we have converted the momentum BSE (with a singular kernel) into a 
homogeneous two-dimensional integral equation for the real weight function p^^[u!,z), i.e.. 



where the kernel V^^^{uj, z; cu, z) is a regular multivariable function. 

The kernel V^^"^ also automatically supports the domain Vt where the function p^^\uj^z) 
is nontrivial. This domain is always smaller then the infinite strip [0, oo) x [—1,1], as is 
explicitly assumed by the boundaries of the integrals over u and z. For instance, with the 
simplest kernel parametrized by a free gluon propagator and constituent quarks of mass m, 
we get for the fiavor-singlet meson p"^\ijj, z) 7^ only ii uj > m?. 

In our approach, to solve the momentum BSE in Minkowski space is equivalent to finding 
a real solution to the real integral equation (12. 6p . No special choice of frame is required. If 
one needs the resulting vertex function, can be obtained by numerical integration over p^v 
in an arbitrary reference frame. 




(2.4) 



F(w, z-P,q) = uj- (g2 + q.Pz + P'^/A) - ie , 



(2.5) 




(2.6) 
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III. NUMERICAL RESULTS 



In this section we discuss the numerical solution of the BSE with various interaction 
kernels. For that purpose, we shall vary the coupling strength as well as the effective gluon 
mass rrig. We are mainly concerned with the range of binding energies that coincide with 
those of heavy quarkonia, which systems we shall study in future work. Moreover, we take 
a discrete set of values for the mass rUg, such that it runs from zero to the value of the 
constituent quark mass. These values are expected to be relevant for the case of a true 
gluon propagator (when nig is replaced by the continuous spectral variable u (12.21) ). Thus, 



in each case, the corresponding gluon density is Pg{c) 
kernel of the BSE to be (in the Landau gauge) 



Ng6{c — m?), which specifies the 



V{q-p) = / 



rrig + le 
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'l' (3.1) 



where the prefactor (including the trace of the color matrices) is simply absorbed in the 
coupling constant. For our actual calculation, we use the bare constituent propagator Si{pi) 
with heavy quark mass M = m (see Appendix A for this approximation). 

Firstly, we follow the standard procedure: after fixing the bound-state mass (v^P^), we 
look for a solution by iterating the BSE for a spectral function with fixed coupling constant 
a = g'^/{47r). Very similarly to the scalar case the choice N = 2 for the power of F in 
the IR of the bound-state vertex function is the preferred one. This choice is a reasonable 
compromise between on the one hand limiting numerical errors and on the other hand 
avoiding the computational obstacles for high N. Here we note that using iV = 1 is rather 
unsatisfactory (comparing with the massive Wick-Cutkosky model), since then we do not 
find any stable solution for a wide class of input parameters g, rrig. In contrast, using the 
value = 2 we obtain stable results for all possible interaction kernels considered here. This 
includes the cases with vanishing rrig, which means that the numerical problems originally 
present in the scalar models are fully overcome here. The details of our numerical 
treatment are given in Appendix B. 

As is more usual in the nonrelativistic case, we fix the coupling constant a = g'^/{An) 
and then look for the bound-state mass spectrum. We find the same results in either case, 
whether P or a is fixed first, noting however that in the latter case the whole integration in 
the kernel K needs to be carried out in each iteration step, which makes the problem more 
computer-time consuming. 

The obtained solutions for varying a and mass rrig, with a fixed fractional binding r] = 
\fP^ I (2M) = 0.95, are given in Table 1. If we fix the gluon mass at rrig = 0.5 and vary the 
fractional binding rj, we obtain the spectrum of Table 2. 



rrig /rrig 


10"3 


0.01 


0.1 


0.5 


a 


0.666 


0.669 


0.745 


1.029 



TABLE 1. Coupling constant as = g'^/lAir) for several choices of rUg/M, with given 

binding fraction r] = v^P^/(2M) = 0.95. 



r] : 


0.8 


0.9 


0.95 


0.99 


a 


1.20 


1.12 


1.03 


0.816 
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0.8 

0.7 



0.6 




FIG. 1: The rescaled weight function r = ^ for the following model parameters: r/ = 0.95, 

rUg = O.OOIM, Us = 0.666; the small mass nig approximates the one-gluon-exchange interaction 
kernel. 

TABLE 2. Coupling = g"^ /{Att) as a function of binding fraction rj = \/P^/[2M), for 

exchanged massive gluon with nig = 0.5M. 

For illustration, the weight function pt^^ is displayed in Fig. [H 

IV. SUMMARY AND CONCLUSIONS 

The main result of the present paper is the development of a technical framework to solve 
the bound-state BSE in Minkowski space. In order to obtain the spectrum, no preferred 
reference frame is needed, and the wave function can be obtained in an arbitrary frame — 
without numerical boosting — by a simple integration of the weight function. 

The treatment is based on the usage of an IR for the Green's functions of a given theory, 
including the bound-state vertices themselves. The method has been explained and checked 
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numerically on the samples of pseudoscalar fermion-antifermion bound states. It was shown 
that the momentum-space BSE can be converted into a real equation for a real weight 
functions p, which is easily solved numerically. The main motivation of the author was to 
develop a practical tool respecting selfconsistency of DSEs and BSEs. Generalizing this 
study to other mesons, such as vectors and scalars, and considering more general flavor 
or isospin structures, with the simultaneous improvement of the approximations (correctly 
dressed gluon propagator, dressed vertices, etc.), will be an essential step towards a fully 
Lorentz-covariant description of a plethora of transitions and form factors in the time-like 
four-momentum region. 
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APPENDIX A: KERNEL FUNCTIONS 

With the Dirac indices explicitly written out, the BSE for quark- ant iquark bound states 
reads 



r(g, PU = ^ I (01^(9 + ^/2)/37r(g, P),,'S{q - P/2)yp,V^f,p,,{q, k- P), (Al) 

where the Lorentz indices of the vertex function have not been specified. 

In our approximation, the IR for a pseudoscalar bound-state vertex function is 

r(g, Pr = if r f dz J'^^'^f (A2) 
^0 J^i [F(w, z- P, q)r 

and the generalized kernel for the BSE in the ladder approximation has the form 



oo 



Vaf3-,5 = 9^J^ dc-^p—^r^f^Pg{c)Yo.pl^^8. (A3) 

where the indices a, /5, 7, 5 and /i, v stand for the appropriate Dirac and Lorentz structures, 
respectively. Moreover, we use the IR for all functions entering the BSE, including the vertex 
[2]), the kernel flA3l) . and the propagators S{q±) (12.11) . 

For the purpose of brevity, we shall use the following abbreviation for the prefactor: 

/•oo r~^^ 

3/ duj dz dc da dbp^^^u, z)g^pg{c) . (A4) 



/<s Jo J-i Jo Jo Jo 
With this convention, the BSE can be written as 



; p'^jCuJ) ^rr dS pMpMiq" - ^V4) - Ps{a)Ps{h) 

[F{CjrzVP.P)f 'JsJ (27r)4 [F(a;,z;g,P)f 
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where the trace is taken over the Dirac indices (after multiplying by 75), and where 

Di = {q + P/2f -a + ie, 
D2 = {q- P/2f -b + ie, 
D3 = {q-pf -c + ie. 



(A6) 



In the following, we shall transform the r.h.s. of the BSE flA5|) into its IR, i.e., the l.h.s. of 
Eq. (|X5|). 

As a first step, we use the algebraic identities 



_uj- q.Pz - P'^/A 
F{uj,z;q,P) ~ F{u,z;q,P) 

g.P _ 1 / 1 1 b-a' 
D1D2 ~ 2[d'i^ D~2^ DiD2^ 

which gives us for the r.h.s. of Eq. (lASp 



(AT) 
(A8) 



N 



+ 



sJ (27r)4 \ [F{uj,z-q,P)f D, \D, D, 
{u - PV2 + -fz)p,{a)p,{b) - Ps{a)Ps{b) 
[Fiu;,z; q, P)f D^D2Ds 



py{a)pr{b) 



[F{uj,z;q,P)f^'DiD2D. 



(A9) 



Furthermore, we employ Feynman parametrization, starting with the first term in expres- 
sion (1A9I) . which yields 



1 



[F{uj,z;q,P)f \D, D 



1 1 



1 , r(N + i)x^"^ 

dx- 







T{N) 



q^ + q.P{zx - (1 - x)) + - uox - a{l - x) 
q^ + q.P{zx + (1 - x)) + - oJX - 6(1 - x) 



~N-l 



-N-l 



(AlO) 



Substituting Eq. (lAlOp back into expression (lA9p . using the Feynman variable y so as to 
match the scalar propagator D^, and integrating over the four-momentum g, we get for the 
first line in expression (1A9P 



\N- 



p2 



1 f pv{a)pv{b) 
Js 2(47r)2 



dx dyy^ x^ ^ z 
Jo Jo 



4 

p2 



y - y\xil + z)- If + _ y)y + Rpy(i _ + z)-l)- U{a) 



y - y\x{z - 1) + \f + q\\ - y)y + P.py(l - y){x{z - 1) + 1) - U{b) 



-N 



where we have defined 



U{a) = [ux + a(l — x)]y + c(l — y) 



(All) 
(A12) 
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Making the substitution x — > 5, such that z = x{l + 2;) — 1 in the second hne of expres- 
sion (lAlip . we can write for the first term in the large bracket ( ]A11[) (including all prefactors 
of the first line (lATTD ) 



1 + z 



z\ N-1 



(4 



where F is defined in Eq. (12. 5p . and where 



N ' 



(A13) 



n 



Introducing the identity 



UJ 



du}5{uo — Vt) 



(A14) 



(A15) 



into the expression ( ]A13ll . changing the integral ordering, and integrating over the Feynman 
variable we arrive at the desired expression for ( 1A13I) 



duj / dz 



r""[F(^,5;p,P)f' 



Is 2(4vr)2 . , . ^.^^ 
where ya± are the roots of the quadratic equation 

y^A + yBa + c = 0, 



with the functions A, Ba, Da defined as 



Of yajd{Da)d{yaj)0{'^ " Vaj] 



(A16) 



(A17) 



A = uj- S] Ba = {uj~ a)T+ + a-c-uj; Da = B^ - AAc; 



(A18) 



Similarly, we can repeat this whole procedure for the second term in expression (lAlip . which 
gives 



T 



N-l 



yb± 



s 2(47r)2 

-Bh ± yp'b 

2A 



-^Pv{a)pv{b)9{z - z)J2 



y,,d{D,)d{y,j)e{l-y,,) 



Bb = {uj- h)T. + b-c-u; Db = B^ - AAc. (A19) 



where we now use the label X2 instead of Xi- 

Next we transform the last term in the second line of expression (lA9p to the desired 
IR form (12 .40 . For that purpose, we basically follow the derivation already presented in 
Ref. which is essentially the same. However, because of the notational differences, we 
give all calculational details below. 
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Let us denote the relevant integral as 



I = i- 



1 



Using the Feynman-parametrization technique, we first write 

^ 1 



(A20) 



dr] 



2_y [M2-/(g,P,r/)-ze]2' 
^„ ; f{q,P,7]) = q^ + r]q.P + Py4. 



a + b a — b 

Then, the denominator of the IR of the bound-state vertex is added: 

1 1 
dr] / dt 



(A21) 



T{N+1) 



[F{uj,z;q,P)]^-^ 2r(iV - 1) 
R = Ljt+ {1 -t)M^ , 



{l-t)t 



N~2 



-1 



[R-f{q,P,z)-te]^+^ 



(A22) 



where z = tz + {1 — t)T]. Now, we match with the function D-^ and integrate over the 
four-momentum q. Thus, we get 



/ = I 



d^q D1D2D3 
(27r)4[F(a.,^;g,Pr-i 

+ j dr^ j dt{l-t)t''-^ j dzx"" I 



(A23) 



2r(iV- 1) 



-1 



d^q 
(2^ 



-q^ + q-Q-{l- x)p^ - + (1 - x)c + xR-ie 



-,~{N+2) 



TiN) 



1 



n{t) 



(47r)2 T{N + 2) - x)^ [n{t) - f{q, P, z) - ie]^ 

R c X 
+ - - 



1 — X X {1 — x) 



S. 



(A24) 



where Q = (1 — x)p — xzP/2 and the function 5* has been defined in Eq. flAlSI) . Note that 
z lies in the interval [—1, +] > and < S* < {^/a + Vb)'^/A. Interchanging the integrals over 
rj and t such that 



1 1 



-10 -1 

T± = and z = tz + {1 — t)r] , 

we finally obtain 



z — z] 



(A25) 



N -I 
2(47r)2 



1 1 

dz 
1 



dx 



(1 -x)^ 



Y.Q{s{z-z)) 



s=± 



dtt 



N~2 



[P(fi(t),5;p,P)] 



N 



(A26) 
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Furthermore, we make the tdependence of F(Q(t), z;p, P) exphcit as 



where 



S- p, P) = ^-p^ t + F(fi(0), ~z- p, P) 
1 — X 



1 — X X 

T^/x -r/ X b + a b — a ^ 
R{t) = J{u,z)t + — —z, 

b + a b — a 
J[UJ, z) = uj — z . 



(A27) 



(A28) 



Integrating over the variable t yields 



/ 



F{m. Pr (N - 1) F(fi(0), ~z;p, P) [F{n{t), z;p, P)] 



(A29) 



and so 



2(47r) 



1 Jo 



dz dx 



9[s{z - ~z)]T, 



N~l 



i {l-x)^F{n{0)J;p,P) [F{A{Ts),~z;p.P)] 



N-l ■ 



(A30) 



In order to separate the F's in the denominator, we use the identity 

1 



F{Q{0),z;p,P) [F{Q{t),~z;p,P)] 



N-l 



1 — X 



1 



1 



J(a;,z)T, [F(fi(0),z;p,P) F(A(T,), 5; p, P) J [F(A(T,), 5; p, P)] 



. (A31) 



Note that, for a given A^, one can repeat this algebra A^ — 1 times till the power of the last 
factor vanishes, which is precisely the reason to use the trick (lA3ip . After this operation, 
the momentum dependence of the denominator in each term becomes formally the same as 
in the desired IR. Although it is possible to derive the corresponding formula for arbitrary 
A^, this would lead to unmanageable expressions (probably not even in closed form), so that 
we rather choose one concrete value of A^. Motivated by the success of the scalar-model 
studies, we take N = 2 henceforth. Explicitly, we obtain for / 



1 



2(47r) 



/ dz dxy^ 
J-i Jo tt 



1 e[s{z~z)] 



J(uj, z){l — x) 



1 



1 



P(fi(0),5;p,P) P(A(T,),5;p,P) 



(A32) 



Integrating by parts over the variable x, we get for the latter expression 



-1 /•! ri „ e[s(z-z)]ln(l-x) 
dz dx y — r 



2(47r)2 7-1 Jo 



dQ{0) 
dx 



dAjTs) 
dx 



F\n{0),~z;p,P) F\AiT,),~z-p,P) 



(A33) 
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Implementing the identity 

1=1 6{Co- Q{t)) (A34) 



into the integrand of expression flA33l) . changing the order of integration, and carrying out 
the integration over the variable x, we arrive at the desired result for the second term in the 
second line of expression flA9p : 



dA{xj{0)) 
dxj{0) 



1 [F(cD,5;p,P)f' 



- Y: Ol'U - S)l '°'-',7'''f'" »|O(0)l»fe(r.)|g|l - x,(T,)]sgn{E{Xi{T.)]}\ , (A35) 

S=± J[UJ,Z) J 

where we have included the previously omitted prefactor, and where the Xj are the roots of 
the delta- function argument in Eq. (]A34p . viz. 



= 2A 

D{T) = B{Tf - AAc , B{T) = R{T) - Cj - c 

^=1^ ,i?(x)=.-5-^. (A36) 

Introducing a more compact notation for the sum over an arbitrary function U of parameter 
T, namely 

^ f/(T) = f/(0) - 9{z - z)U{T+) - e{z - z)U{T^), (A37) 
r 

we can rewrite Xs in Eq. flA35l) as 

UCo,~z) = I E E - x,{T)HnE[x,{T)]. 

(A38) 

The first term in the second line of expression (1A9P can be treated in a very similar fashion 
as the previous case, though accounting for the different power of F in the denominator. 
Doing so explicitly, and multiplying with the correct prefactor, the resulting expression reads 



(A39) 

Assuming validity of the Uniqueness Theorem, we have converted the momentum BSE into 
an equation for the weight function. It reads 



p^'^\Co,z)= duo dzV^'^\u,z]uo,z)p^^\u,z) , (A40) 
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where the kernel is simply given by the sum of the contributions derived above: 

/•oo /-oo /'OO ^ 

V'-'^'{uj, z;u!, z) = —Sg"^ / dc da (i6p(,(c) Xj(u), 5; cu, 2;). (A41) 

Jo Jo Jo ~[ 

1. Heavy quark approximation — unequal- mass case 

When the quark is sufficiently heavy (say mq{2 GeV) >> Aqcd)) then the approximation 
of the quark-mass function by a constant turns out to be adequate. Neglecting self-energy 
corrections is equivalent to the use of a free heavy-quark propagator, which corresponds to 
the free-particle spectral functions 

Mip^(a) = ps{a) = M,5{a - M^) , 

M2p,{b) = ps{b) = M^Sib - Mi) , (A42) 

where the variables a, b distinguish the types of quarks the bound state is composed of. For 
completeness, we write the kernel down explicitly for this case: 

y P] z; z) = dcpgic) ixi + X2 + X3 + X4) ; (A43) 



1 + 2 ± il-ya±)VDa 

_ ^ zOjz - z) ^ yb±e{D,)e{y,±)e{l - y,±) 



1-2 ± {l-yb±)^/Db 

X3 = E E ^''^K'i ^jr^ d[D{Tmx,{T)]e[l - x,(T)]sgnE[x,(T)], 
r j=± Jv^.z) 

( P2 Ml -Ml \ 

g[Z^(T)]g[x,(T)]g[l-x,(T)] r rj(^,2) 1 _ , .^^.l 

^,ti ^^(^,2) l(l-x,(r))|E[x,(r)]| x,(7))s.ni.[x,(i)J|. 

(A44) 

Here, the arguments a, 6 of the functions x, J, D must be replaced by the quark masses 
Mi,M2, respectively. 



2. Equal-mass case 



In the case of quarkonia, the kernel becomes more symmetric with respect to the variable 
2, so that the formula for the kernel further simplifies. The function R depends on z only 
through the variable T, such that 

J{uj,z) J = uj- M^; R{T) ^ R{T) = JT + M^, (A45) 
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where M is the common mass M = Mi = M2. The roots then become 

ya± Xj=±(T+); yb± Xj=±{T_) , (A46) 

and the kernel can be written in a more compact form: 

V^^\u,z;u;,z) = H dcp,{c)K^'\u, z;u;, z, c); (A47) 

2(47rj^ Jo 

K'^^ (a), z] u, z, c) = 



s=±j=± (1 + sz)^D{Ts){l - Xj{Ts)) 

T j=± 

~ ' ln(l - x,(T))sgnE[x,(T)]^ \ . (A48) 



3. One-gluon-exchange approximation 

In order to consider the one-gluon-exchange approximation, we must take the massless 
limit c — > and restrict ourselves to the equal-mass case. One can easily recognize that the 
root x+ = becomes trivial and the associated contribution in expression (1A48P vanishes. 

For the purpose of brevity, we label 

x(r)^x_(T) = ^:^f^ (A49) 
Taking into account the relations 

v^/x = A, E = A, 
and doing a little algebra, we get for the kernel 



Y^e[s{z-~z)]e{-Bs)e{A + B, 

s=± 



A + B, V A / 4J2 



+ In 1 + 



(A50) 

For completeness, we recapitulate here the complete list of functions: 

A = Cj-S, J = lj-M\ S = {1-z'^) — , 
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4. Case N = 1 



Repeating the derivation, but now for the parameter value = 1, we should obtain the 
homogeneous equation 

pW(cD,5) = y^ du J ^dzV^^\Cd,z;Lj,z)p^^\u,z), (A52) 

where the kernel is given by the expression 



V^^\Cu,z;uj,z) 



s 2(47r) 



2 



X2 = Pv{a)p^{b)- e{z -z)}_^ 

i- — Z V 

X.-P.[a)p^ib)^^ -^^^ , 



_ d[D{T)]e[x,{Tmi-x,{T)] 
Xa- l^ — 



\E[x,{T)]\J{uj,z) 



(A53) 



Here we use the same notations and conventions as in the previous section. The correspond- 
ing derivation is very straightforward and exactly repeats the steps made for the case N = 2, 
so we only add a few comments. The expressions for xi.2 have in fact been derived in the 
previous section, as they were for general A^ see Rel. (lA16p and Rel. (lA19p . The expression 
for Xa we adopt from Ref. while the remaining function Xs follows from conversion of 
the term with F° = 1, i.e., 

' (A54) 



which is present in the corresponding momentum-space BSE. The proper derivation is in 
fact easier than in the N = 2 case, and the result represents the basic PTIR for a scalar 
triangle with one leg momentum constrained such that < {Mi + M2Y. 



APPENDIX B: NUMERICAL PROCEDURE 



In this appendix we describe the numerical procedure actually used for obtaining the 
bound-state spectra. Because of the structure of the integral equations to be solved, the 
treatment is similar to the procedure used in Refs. [l0|,[li|. However, we do introduce a rather 



tricky modification here, which improves the numerical stability when the mass parameter 
mg is too small as compared to the constituent mass. Details of this technical difference are 
given below. But first we describe the numerical treatment for the case mg ~ M [M is a 
heavy quark mass). 
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Equation (1A48P is a homogeneous linear integral equation whose solution needs to be 
properly normalized. For that purpose, let us adopt the auxiliary normalization condition 

i=rdzrdj-^^. (Bi) 

J-i Jo 

Then we find that the equal-constituent-mass BSE can be transformed into the inhomoge- 
neous integral equation 

(Co, z) = Kf {Co, z) + duo J^^ dzK'§ {Co, z; uo, z)p^^^ {uj, z) (B2) 

with the kernels 

kP{u,S) = ^^[D(0)]^K(0)]^[1 -a;,(0)]ln(l -x,(0))sgnE[x,(0)], (B3) 



K^^\^,~z;u;,z,c) = ^ E g[D(0)]g[a;,(0)]g[l - x,(0)] ln(l - x,(0)) 'g^^^^^j°^J ^ 



+ ^ E EjH^ - mD{Tsmx,{zmi - x,{t,)] x 

sz TsXj{Ts) ^ 



l + sz)JD{Ts){l-x,{T,)) 



El 

2J 



^ ln(l ^^■(^-)) sgnE[x,(T.)]-^. (B4) 



|E[x,-(r,)]|(i-x,(r,)) J ^ I ^'^2j 

The definitions of all these functions have been given in the previous appendix; for the 
readers' convenience we recapitulate them here: 



-B{T) ± JD{T) 
x^{T) = ; D{T) = B{Tf - AAm, ; 



a ' 



A = Co-S- B{T) = R{T) -Co-m 

E{x) =Co I--S; R{T) = JT + ; J = uo - M'^ ; 

x^ 

S=(l-?):^iT, = i||. (B5) 

Equation (lB2p has actually been used for the numerical solution of the case rrig/M ~ 1. 

The kernel in Eq. (IB2p is free from running singularities owing to the presence of G 
functions. However, in the case of a massless-gluon kernel, we would have a singularity just 
on top of the boundary. There, J — > as approaches the quark mass. We find that this 
instability is avoided if we generate the inhomogeneous term in the following manner. First, 
we add exactly zero, in the form 

n^'^Uu z) 

{f{Co,~z)-f{Co,z))P-^^, (B6) 
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to the r.h.s. of the original homogeneous BSE for the weight function p. Then, solving the 
equation 



doj I dz 



with the normalization condition 

1 



1 







(B7) 



(B8) 



is equivalent to the solution of the original BSE. As a suitable function, we choose one with 
the property 

f{u,~z) = l- (B9) 

We observe that this method is applicable for any positive value of m^, but is only slowly 
convergent when it is used for the previously discussed case nig ~ M. We also note that, up 
to a small numerical error, Eq. flB2p yields the same spectra for those a's and m^'s which 



allow both Eqs. flB2IIB7p to be numerically stable. 

Equation flB7l) is solved by the method of iteration. When straightforward iteration 
fails — measure being a difference between weight functions obtained in different iteration 
steps, and/or a deviation of the auxiliary normalization integral from a predefined value — 
we change the coupling constant (in the treatment with fixed P^, otherwise the procedure 
is the converse) until a solution is found. For the numerical solution we discretize the 
integration variables u and z using Gauss-Legendre quadrature, with a suitable mapping 
(— 1, +1) (cUmin, oo) for u . Equatious (]B7|B2p are solved on a grid of = A''^ * A"^ points 
spread over the rectangle (—1, +1) x (cUmin, oo). The value Umin is given by the support of 
the spectral function. In all cases we take A^ = 2A^. Examples of numerical convergence 
for some bound-state cases are presented in Table 3. As we can see, there is a rather weak 
dependence of the eigenvalue a on the number of mesh points A^. The last column gives 
values calculated from a weighted average (WA), with A^^^ the appropriate weight. 





32 


40 


64 


80 


WA 


rj = 0.95; mg/M = lO"^ 


0.6611 


0.6690 


0.6697 


0.6734 


0.669 


r] = 0.95; nig/M = 0.5 


1.037 


1.0259 


1.0210 


1.0229 


1.029 


r] = 0.99; nig/M = 0.5 


0.818 


0.8127 


0.8158 


0.8155 


0.816 



TABLE 3. The coupling = (7^/(47r) for the ladder BSE with fixed ratio rric/M, as a 

function of the number of mesh points. 
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